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NONPARAMETRIC  BAYESIAN  BIOASSAY 


INCLUDING  ORDERED  POLYTOMOUS  RESPONSE 

by 

Alan  E.  Gelfand  Lynn  Kuo 

Department  of  Statistics,  U-120 
University  of  Connecticut 
Storxs,  Connecticut  06269—3120,  U.S.A. 

SUMMARY 

Previous  attempts  at  implementing  fully  Bayesian  nonparametric  bioassay 
have  enjoyed  limited  success  due  to  computational  difficulties.  We  show  here 
how  tnis  problem  may  be  generally  handled  using  a  sampling  based  approach  to 
develop  desired  marginal  posterior  distributions  and  their  features.  A  useful 
extension  is  presented  which  treats  the  case  of  ordered  polytomous  response. 
Illustrative  examples  are  provided. 


Key  Words:  bioassay,  Dirichlet  process  prior,  Gibbs  sampler,  potency  curve, 
product  Beta  prior. 
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1.  Introduction 

We  first  formalize  the  quantal  bioassay  problem  as  follows.  An 
experimenter  wishes  to  investigate  the  potency  of  a  stimulus  by  administering  it 
at  k  dosage  levels,  tp  ^  tota*  °*  ni  su^iects  are  treate^  a*  level  t 

with  the  number  of  positive  responses  obtained  denoted  by  Xj,  i=l,2,...,k.  The 
potency  curve  F  is  the  distribution  of  tolerance  levels,  that  is  F(t)  is  the 
probability  of  achieving  a  positive  response  at  dosage  level  t.  We  seek  to  make 
inicreuu;  about  F.  Parametric  models  as,  for  example,  discussed  in  Finney 
(1978)  typically  specify  F  as  a  family  of  logit  or  probit  forms.  We  take  a 
nonparametric  setting  only  assuming  F  to  be  an  arbitrary  right-continuous 
non-decreasing  function  whose  range  is  [0,1]. 

Taking  a  Bayesian  inferential  framework,  we  note  that  the  fully  Bayesian 
nonparametric  approach  to  estimating  the  tolerance  distribution  in  a  quantal 
bioassay  has  not  previously  been  implemented  due  to  computational  difficulties. 
Antoniak  (1974)  has  shown  that,  under  Ferguson’s  Dirichlet  process  prior  (1973), 
the  posterior  distribution  of  the  potency  curve  is  a  mixture  of  Dirichlet  process 
distributions.  Unfortunately  this  mixture  becomes  increasingly  intractable  as  the 
number  of  observed  dosage  levels  increases.  In  her  unpublished  dissertation,  of 
M.N.  Wesley  abandons  calculation  of  the  posterior  expectation  of  the  potency 
curve  at  the  observed  dosage  levels  because  of  such  difficulties,  estimating, 
instead,  the  posterior  mode.  Kuo  (1988),  again  in  deference  to  these  difficulties, 
obtains  linear  Bayes  estimates  under  squared  error  loss.  Related  previous 
Bayesian  work  includes  Kraft  and  Van  Eeden  (1964),  Rimsey  (1972), 
Bhattacharya  (1981),  Disch  (1981)  and  Ammann  (1984). 
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This  paper  has  two  objectives.  First,  we  show  how  fully  Bayesian  analysis 
to  obtain  for  any  t,  the  posterior  distribution  of  F(t)  given  the  data  X  =  (Xp 
X2,...,Xk)  can  be  straightforwardly  implemented  under  two  rich  classes  of  prior 
specification.  In  addition,  features  of  these  distributions  such  as  mode, 
expectations  and  quantiles  can  be  readily  obtained.  Second,  we  provide  a  useful 
extension  of  the  quantal  bioassay  model  which  allows  ordered  polytomous 
response  arising  from  stochastically  ordered  potency  curves.  Missing  data  can  be 
readily  handled.  To  our  knowledge  no  literature  on  Bayesian  approaches  to  this 
problem  exists. 

The  format  is  the  following.  In  Section  2  we  develop  the  details  to 
implement  the  first  objective.  Section  3  provides  an  illustrative  example  using 
the  Dirichlet  process  prior  with  a  data  set  from  Cox  and  Snell  (1989).  In  Section 
4  we  develop  the  aforementioned  extension  while  in  Section  5  we  provide  an 
illustrative  example  adapted  from  Maxwell  (1961)  for  the  case  of  two  ordered 
response  levels  utilizing  a  product— Beta  class  of  priors.  Section  6  offers  brief 
summary  and  discussion. 

2.  Models,  Distribution  Theory  and  Inference  For  the  Basic  Bioassay  Problem 

We  assume  the  responses,  X-,  to  the  dosage  levels  t-  are  independently 
distributed  as  Binomial  Bi(n- ,p- )  with  p^  =  F(t-)  where  F(t)  is  an  unknown 
underlying  potency  curve  yielding  the  probability  of  a  response  at  dosage  level  t. 
Hence  the  likelihood  at  X  =  x  is 


x 


L(p;x)  =  n  (XJ)  pj1  (l-pj) 


n  — x 

i  i 


«  A  ■ 
1  =  1  1 


(1) 
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where  p  =  (p1,...,pk). 

Interest  focuses  on  inference  regarding  the  p.’s  and  the  function  F. 

Often  F  is  given  a  parametric  form  F(t;0)  such  as  a  scale  and  location  logistic 
or  probit  curve.  This  literature  is  large.  See,  for  example,  Finney  (1978)  for 
discussion  and  references.  Here  we  assume  F  belongs  to  the  nonparametric  class 
of  right-continuous  non— decreasing  functions  taking  values  in  [0,1]. 

The  Bayesian  framework  requires  specification  of  a  prior  which  represents 
our  beliefs  regarding  F(t).  Practically  this  necessitates  specification  of  a 
joint  measure  for  p  with  density  denoted  by  h(p)  wuicn  by  assumptions  on  F 
is  over  the  set  S  =  {p:  0<pj<...<pk<l}.  One  family  of  priors  which  has  been 
discussed  in  this  context  in  the  literature  is  the  Dirichlet  process  prior  introduced 
by  Ferguson  (1973).  Kuo  (1988)  provides  a  recent  summary  of  this  discussion. 
This  specification  assumes  that  F  is  close  to  some  given  FQ  with  closeness 
quantified  by  a  given  precision  M.  More  precisely,  for  any  t,  the  induced  prior 
on  F(t)  is  a  Beta  distribution,  Be[MFQ(t),  M{l-F0(t)}].  Since  E{F(t)}  = 

FQ(t)  and  Var{F(t)}  =  F0(t){l-F0(t)}/(M+1),  Fq  and  M  have 
interpretations  which  facilitate  their  specification  from  prior  information.  In 
particular  Fq  is  usually  taken  to  be  a  standard  distribution  whose  median  agrees 
with  our  prior  guess  for  the  LD50  and  whose  spread  provides  rough  agreement 
with  our  prior  expectations  at  other  dosage  levels.  M  is  chosen  to  reflect  our 
confidence  in  Fq  in  accord  with  the  extent  of  our  prior  experience.  In  practice 
we  might  experiment  with  several  choices  of  M  to  see  the  effect  on  posterior 
features  of  interest. 

The  induced  prior  on  p  is  an  ordered  Dirichlet 
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k+1 

'Dfri-k+r1 1  -pi71  (P2-P1)72  *•  (pk-Pk-i)7k  1<i-pk)7k+i  1 
n  r(r) 

i  =1  1 


(2) 


where  7^  =  M{FQ(t.)  -  i=l,...,k+l.  Here  we  assume  tQ  =  -®,  F  (t  ) 

=  lk+l  =  0D’  ^o^k+l^  =  so  ^i  =  Marginal  and  conditional 
distributions  and  expectations  under  (2)  are  routine. 

A  second  class  of  priors  for  p  is  the  product— Beta.  This  family  has  been 
discussed  in  unpublished  work  of  L.  Sharpies  and  takes  the  form 

k  a—1  0.-1 

*B(l)=ckM  n  p.  (1— Pj)  ,  oi>0, /?j>0  (3) 

where  a  =  (a1,...,Qk),  0  =  {0^,...,0^)  and  ck  is  the  normalizing  constant  under 

k 

restriction  to  S  .  Sharpies  shows  that  ck  can  be  expressed  as  a  finite 
multidimensional  summation. 

Note  that  Xg  offers  mathematical  convenience  in  that  it  is  conjugate 
with  respect  to  (1)  while  Xg  is  not.  The  family  (3)  is  a  flexible  class  of  priors 
but  it  is  not  clear  how  to  select  a  and  0  in  accord  with  prior  information. 

That  is,  unlike  Xg,  Xg  is  not  induced  in  any  obvious  way  from  a  prior  on  F 
since  it  is  specified  within  dosage  levels  rather  than  across  them.  However,  Xg 
can  be  chosen  to  reflect  to  prior  information  about  F  in  terms  of  a  given  Fq 
and  given  precision.  Letting  e^  denote  a  row  vector  having  value  1  at  the  i^ 
coordinate,  0’s  elsewhere,  expectations  involving  the  p.  may  be  formally  given, 
e  g.,  EsCPj)  =  ck(Q>  0)/ck(o  +  e^,  0).  This  suggests  equating  Efi(pi)  = 

F0(tj),  i  =  l,...,k.  Moreover,  Mj  =  Oj  +  0^  can  be  viewed  as  a  precision 
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parameter  analogous  to  M  above.  The  magnitude  of  NT  reflects  our  confidence 
in  the  value  FQ(tj).  With  specification  of  hi  we  obtain  k  equations  in  k 
unknowns.  Unfortunately,  explicit  calculation  of  the  Eg(p.)  will  be  infeasible 
except  in  very  special  cases  making  solution  of  this  system  of  equations  virtually 
impossible.  However,  the  conditional  distribution  of  p.  given  pj,  j#  is 

obviously  a  Beta  distribution,  Be(Oj, /3j)  restricted  to  [p. _ Pj+1]-  Though, 

again,  the  mean  of  this  conditional  distribution  will  not  be  available  explicitly 
the  mode  is  readily  obtained  It  is  p j  =  (a— 1)/(M.— 2)  provided  a->l,  /?.>1 
and  p-_j  <  p.  <  Pj+j-  i axmg  />•  as  an  approximation  to  the  marginal  mode  for 
p.  we  equate  it  to  F0(t)  whence  a  =  (Mj-2)  F0(t)  +  1  and  /?•  =  hi  -  a.. 

Objects  of  primary  interest  for  Bayesian  inference  are  the  marginal 
posterior  distributions  of  Pj|X  and,  at  a  specified  t,  of  F(t)|X.  For  such 
distributions  the  mean  or  the  mode  provide  point  estimates  while  appropriate 
quantiles  provide  interval  estimates.  As  noted  in  the  introduction,  computation 
of  such  distributions  and  estimates  has  proved  very  difficult.  However  the 
sampling  based  approach  discussed  in  the  context  of  hierarchical  Bayes  models  in 
recent  papers  by  Gelfand  and  Smith  (1990),  and  Gelfand  et.  al.  (1990)  is  ideally 
suited  to  this  problem.  This  approach,  known  as  the  Gibbs  sampler,  is  an 
iterative  Markovian  updating  scheme  which  dates  at  least  to  Metropolis  et  al 
(1953).  We  do  not  review  details  here  merely  remarking  that  implementation 
requires  sampling  from  so  called  complete  conditional  distributions.  For  the 
remainder  of  this  section  we  develop  these  distributions  under  both  Xg  and  x^ 
and  indicate  how  the  desired  posterior  density  estimates  and  features  are 
obtained. 

We  note  that  in  our  illustrations  sampling  is  conducted  with  v 
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independent  parallel  replications  each  taken  to  r  iterations.  Choice  of  v 
determines  how  close  our  density  estimate  is  to  the  exact  density  at  the  r1*1 
iteration  the  order  of  convergence  being  0(v— *).  Choice  of  r  determines  how 
close  the  latter  density  is  to  the  actual  marginal  posterior  density  with 
convergence  at  an  exponential  rate  (Geman  and  Geman,  1984;  Tanner  and  Wong, 
1987).  Settings  for  r  and  v  to  achieve  smooth  converged  estimates  vary  with 
the  application  and  require  diagnostic  assessment  as  in  Gelfand  et.  al.,  (1990). 

For  the  examples  in  sections  3  and  5  r  =  20  and  v  =  1000. 

Each  of  the  priors  (2)  and  (3)  requires  a  minor  ivust  to  facilitate 
implementation  of  the  sampling  approach.  For  ir^,  although  (1)  and  (2)  are  not 
conjugate  with  respect  to  the  p-,  introduction  of  a  set  of  unobserved  multinomial 
variables  simplifies  the  required  sampling.  Let  Zj  =  (Zjp...,  Zi,k+l)  • 

Mult(n-,  A),  i=l,...,k,  where  A=(Alf...,Ak+1)  with  A j  =  Pj  “  Pj_i>  P0=  °>  Pk+1 
=  1.  The  variable  Z-j  denotes,  amongst  the  n-  individuals  receiving  dosage 
level  t.,  the  unobserved  number  who  would  have  responded  to  dosage  level  tj 
but  not  to  dosage  level  t^_j. 

The  Gibbs  sampler  may  be  implemented  using  random  draws  from  the 
complete  conditional  distribution  of  p|X,  Zj,...,Zk  and  from  the  complete 
conditional  distributions  of  Z- 1 X,  p,  Zj,  jft.  The  former  is  an  ordered  Dirichlet 
updating  (2) 


k+1 


^A  V  k+1  Tj— 1 

n  r(7  )  J-1 

j=i  J 
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where  7j  =  7j  +  £  Zy.  Thus,  the  complete  conditional  density  for  p. ,  over  the 
set  (Pj_i»  Pi+1]>  denoted  by  gD(Pj  I X,  Z1,...)Z]£,  pj,  tfi),  is  that  of  ^  +  pf_j 

where  -  (p^^Pj.^BetaC  7^  ^i+1)- 

The  complete  conditional  distribution  for  Z.  is  a  product  of  two 
multinomials.  That  is,  writing  Z-  =  (Zj(l),  Z-(2))  where  Z-(l)  =  (Z^,...^.), 
Zj(2)  =  (Zj  j+1,...,Zj  Z^l)  and  Zj(2)  are  conditionally  independent  given 

Xj  with  Z.(l)  -  Mult(X.,  Pj1  A(l)),  Zj(2)  -  Mult(nrX.,  (l-pi)"1A(2)),  A(l)  = 
(Ap...,Aj),  A(2)  =  ( A- ^ 2 )- 

After  v  independent  replications  each  to  the  r^  iteration  we  obtain 
(p(r),  Z^,...,zj^),  s=1v»v-  The  resulting  estimate  for  the  marginal  posterior 
density  of  pj  is 

WX>  =  V-1  j^otPjix,  pjj),  j#i)  (4) 

Similarly  the  posterior  mean  of  p^  is  estimated  using  the  mean  of  g^  leading  to 


ED(Piix)  =  »  IszipW)5+(p(;],rp!:),sH7is/(Tisn+1,s)})  <5> 


where  7JS  =  7j  +  £  z|j]»  j=l,...,k+l. 

A  posterior  density  estimate  for  F(t)  at  say  t=t  may  be  obtained  by 

♦ 

including  F(t  )  as  an  additional  model  parameter.  More  precisely  we  revise  the 

*  * 

prior  to  include  F(t  )  and  to  take  the  form  ^(p)  •  h^(F(t  )[p).  Paralleling 

*  *  * 

(2)  if  t  t  [t- ,  tj  jJ,  h  is  naturally  the  density  of  p;  +  a  where  a  -  (Pj+1  - 
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*  *  *  * 

Pj)  Beta(7  ,  7i+1  -  7  )  with  7  =  M{FQ(t  )-F'0(ti)}.  Since  there  is  no  data 

*  * 
at  dosage  level  t  ,  the  complete  conditional  distribution  for  F(t  )  is 

*  * 
hjJF(t  )  j  p)  as  well.  Therefore  the  posterior  density  estimate  for  F(t  )  is 

—1  *  M 

v  E  hp(F(t  )|p^  ').  A  convenient  by-product  of  this  specification  is  that  the 

* 

posterior  mean  of  F(t  }>* 

ED{F(t*)|X}  =  Ep  {Pj  +  (pj+1  -Pj)  7*/7i+1|X} 


* 


Expression  (6)  shows  that  once  the  posterior  means  for  the  p- *s  have  been 
computed  using  (5),  an  elementary  interpolation  formula  enables  Eg(F(t)  |  X) 
for  any  t.  This  formula  was  first  noted  by  Antoniak  (1974).  The  unpublished 
Ph.D.  theses  of  M.N.  Wesley  also  discusses  (6)  as  well  as  computation  of 
Ep(p-  |X)  In  fact,  she  obtains  exact  formulas  for  these  eructations  for  UP  to 
four  dosage  levels  but  evaluation  of  these  formulas  even  for  three  levels  is  a 
forbidding  task.  By  comparison,  the  estimate  (5)  is  routine  to  evaluate 
regardless  of  the  number  of  dosage  levels.  Moreover,  using  the  above  estimates  of 
the  posterior  densities  for  the  p^  and  F(t)  we  may  straightforwardly  obtain 
other  features  of  these  densities  such  as  modes  and  quantiles. 

For  rg,  examination  of  (1)  and  (3)  reveals  that  the  complete  conditional 
distribution  of  p^  |  X,  pj,  jri,  is  that  of  a  Beta(a  +  X^  ^  ^  -  Xj)  restricted 

t0  [Pj_i-  Pj+il  which  we  denote  by  ggtpJX^  p^j,  pj  +  1).  Sampling  from  gB 
by  retaining  appropriately  restricted  draws  from  the  unrestricted  Beta 
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distribution  will  be  very  inefficient  when  —  p-_^  is  small.  One— for-one 
sampling  from  gg  may  be  achieved  through  the  use  of  the  incomplete  Beta 
mnction  which  is  incorporated  into  many  scientific  subroutine  libraries.  If  G(y) 
=  P{Y<y|  Y-Beta(a.b)}  then  Z  is  distributed  as  Beta(a,b)  restricted  to  [c,d]  if 
Z  =  G-1[G(c)  +  U{G(d)  -  G(c)}]  where  U  is  a  U(0,1)  variate.  After  v 
replications,  each  to  the  r^  iteration  we  obtain  p[r),  s=l,...,v.  To  obtain  a 
posterior  density  estimate  for  p-  we  would  compute 


fB(PilX)  =  V  1  pS,s.  Pili,s) 


(?) 


Use  of  (7)  is  preferred  to  a  kernel  density  estimate  based  upon  the  p(^  since,  by 
using  the  conditional  structure  of  the  model,  substantially  smaller  v  is  needed 
(see  Gelfand  and  Smith,  1990).  Though  standardization  is  required  to  obtain 
each  gg  in  (7)  this  requires  only  a  univariate  numerical  integration  and  can  be 
done  quite  rapidly  using  simple  trapezoidal  or  Simpson's  rule  integration.  Since 
moments  of  gg  are  not  available  explicitly,  an  expression  such  as  (5)  is 
unavailable  for  Eg(p- 1 X);  posterior  moments  of  Pj  are  most  easily  calculated 
using  the  p[^. 

* 

A  posterior  density  estimate  for  F(t)  at  say  t=t  may  be  obtained  by 

* 

including  F(t  )  as  an  additional  model  parameter.  More  precisely,  as  with  t 

*  .  * 
we  revise  the  prior  to  include  F(t  )  and  to  take  the  form  flg(p)  •  hg(F(t  )|p). 

* 

Paralleling  the  development  below  (3)  if  t  e[t- ,  t  +  1],  hg  is  naturally  a 

*  *  *  *  *  * 

Beta(a  ,0  )  restricted  to  lPj>Pj+i)  wherewith  M  =  a  +  0  specified,  a 

*  *  *  #  * 

=  (M  -2)FQ(t  )  +  1  and  0  =  M  —  a  .  Again,  in  the  abponce  of  data  at 
*  *  * 
dosage  level  t  ,  the  complete  conditional  distribution  for  F(t  )  is  hB(F(t  )|p) 


11 


so  that  the  posterior  density  estimate  for  F(t  )  becomes  v-1  E  hg(F(t  )  |  PgF^) 
Again  one  dimensional  numerical  integration  to  standardize  hg  is  required. 
Unfortunately  no  interpolation  formula  analogous  to  (6)  is  available  for 
EB(F(t)|X). 

3.  A  Numerical  Example 

The  data  for  this  example  are  taken  from  an  illustrative  table  in  Cox  and 
Snell  (1989  p.7).  In  all,  there  are  150  subjects  at  5  different  stimulus  o0 
subjects  at  each  level.  The  data  are  given  in  Table  1 .  A  probit  model  N(2,.25) 
was  chosen  as  the  prior  shape  Fq  with  four  choices  for  strength  of  prior  belief, 

M  =  1,  30,  100,  300.  Table  1  supplies  the  maximum  likelihood  estimate  and 
prior  guess,  F0(t;),  for  each  stimulus  level.  Using  the  Xg  prior,  (2),  Table  1 
also  shows  the  nonparametric  Bayes  estimates,  (5),  of  F(t)  for  each  value  of  M 
using  the  Gibbs  sampler  with  a  total  of  r  =  20  iterations  and  v  =  1000 
replications.  The  supplied  standard  deviation  (S.D.)  measures  the  variation 
amongst  these  replications  and  also  is  an  estimate  of  the  marginal  posterior 
standard  deviation  of  F(t).  Interval  estimates  for  the  F(t- )  were  obtained 
using  the  empirical  distribution  of  the  1000  estimates  from  the  1000  replications. 
In  particular,  equal  tail  95  percent  confidence  intervals  are  given  here. 

Under  Xg  the  Bayes  estimate  of  F(t)  can  be  obtained  using  the 
interpolation  formula  (6).  Figure  1  graphs  the  Bayes  estimate  of  the  entire 
potency  curve  using  this  formula  for  M=1  and  M= 100  as  well  as  the  curve  F 
With  increasing  M,  the  Bayes  estimate  approaches  the  prior  F  .  Figure  2  plots, 
for  M  =  1,  the  Bayes  estimate  with  95  percent  posterior  confidence  band 
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developed  from  pointwise  interval  estimates  using  five  points  between  each 
observed  stimulus  level.  With  increasing  M  the  bands  would  become  narrower. 

4.  Ordered  Potency  Curves 

We  extend  the  model  of  Section  2  to  allow  ordered  polytomous  response 
arising  from  stochastically  ordered  potency  curves.  Suppose  for  a  given  stimulus 
we  observe  whpther  or  not  each  of  two  events  occurred  such  that  one  is  contained 
within  tuc  oiuei.  For  example,  event  A  might  be  "patient  died"  while  event  C 
might  be  "patient’s  condition  worsened"  whence  A  C  C.  Other  illustrations 
include:  in  response  to  a  dosage  of  medication,  A  is  "improvement  in  a  particular 
metabolic  index",  C  is  "improvement  or  no  change  in  this  index";  in  response  to 
hours  of  tutoring,  A  is  a  "high  pass"  on  a  standardized  exam,  C  is  a  "pass"  on 
this  exam.  We  show  here  how  the  approach  of  Section  2  can  be  extended  to 
handle  such  situations  with  further  extension  to  more  than  two  nested  events 
being  obvious.  The  case  of  non-nested  events  as  for  example  if  A  is  a  decline  in 
blood  pressure  while  C  is  a  decline  in  cholesterol  level  is  not  handled  here. 

Formalizing  notation,  suppose  at  level  t- ,  i=l,2,...,k  we  observe  m 
subjects  with  event  A  occurring  Xj  times,  event  C  occurring  Y-  times  with  A 
C  C  so  that  X-  <  Y-.  We  model  this  situation  with  two  underlying  potency 
curves  F^(t),  F^(t)  which  are  stochastically  ordered,  i.e.,  F^(t)  <  F^,(t). 
Letting  Pj=  FA(ti)  <  qj  =  F^ftj)  we  assume  that  the  joint  distribution  of  X; 
and  Yj  is  specified  through  (Xj,  Yj— Xi )  -  Mult  (n^  pj(  qj-p.,  1-qj).  Hence 
with  p  =  (Pp-.-.p^)  and  q  =  (q^.-.q^)  the  likelihood  at  X  =  x,  Y  =  y  is 
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k 

n 

i=l 


n- ! 
1 


xjT(y^r]T(nj^)I 


(8) 


Interest  focuses  on  inference  regarding  the  p^s  and  q^s  as  well  as  and 
F^-,.  Prior  specification  is  restricted  to  the  set  T  =  {(p,q)  :  0  <  Pj  <...<  p^  <  1, 
0  <  qj  <  q2  <  <  q^  <  I,  p-  <  q-  for  all  i}.  We  consider  two  families  of  priors 
extending  and  jrg  respectively. 

In  extending  the  Dirichlet  process  prior  is  replaced  by  a  product 
Dirichlet  process  prior  with  stochastic  :>rdei\  We  do  not  attempt  formal 
definition  or  investigation  of  such  measures  here.  Rather,  if  F^  is  close  to  some 
given  F^  Q  with  precision  and  F^  is  close  to  some  given  F^,  Q  with 
precision  then  for  any  t  we  assume  that  the  joint  prior  on  (F^(t),  F^,(t)) 
is  of  the  form 

M .  F .  f t)-l  M  .  { 1-F .  rt(t)}-l 

4(‘){FA(t)}  A  A’°  {l-FA(t)}  A  A’° 


MrFr  n(t)-l  Mr{l-Fr  n(t)}-l 

{Fc(t)}  C  C’°  {1— Fc(t)}  C  C’° 


(9) 


.  over  0  <  F^(t)  <  F^t)  <  1.  Expression  (9)  is  a  product  of  Beta  forms 
standardized  by  d(t)  under  the  order  restriction.  Similarly  the  induced  prior  on 
(p,  q)  over  T  takes  the  form  of  a  product  of  ordered  Dirichlets, 

k+1  7.-I  77.-I 

nD(p,q)  =  c(7, v)  n  a.j  e.J 

v  *  *  j=i  J  J 

where  Aj  =  Pj  -  Pj_j,  <J  =  qj  -  qH.  Tj  =  MAfFA,o(lj)  -  FA1o<,H))’  Vj  = 


(10) 
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Vco^^Co^H11  an<*  is  the  standardizing  constant. 

7Tg  is  extended  to  produce  a  form  over  T  which  is  conjugate  with  (8), 

that  is, 


*fi(P,q)  =  c(a,  0,  6)  n  p. 
B  j=i J 


Qj  1 


(qj-Pj)  J 


/3;-l  6.- 1 

(1-*)  J 


J 


(11) 


where  again  c  is  the  standardizing  constant.  Linking  5Tg(p,q)  to  prior 
information  can  be  done  similar  to  that  for  flg(p)  using  the  condition?']  modes 
of  Pj  and  of  qj-Pj-  Assuming  a-  +  0-  +  &  =  Mj,  Mj  specified,  with  a.,  /?.,  6- 
>  1  we  obtain  Oj  =  (Mp3)  F^t),  (3-  =  (Mp3)  (F^ri)  -  FA  0(tj)}  and 


=  Mi-ai-/3i. 

Implementation  of  our  sampling  approach  for  the  prior  (10)  is  simplified 

by  the  introduction  of  unobserved  multinomial  variables  as  in  Section  2. 

A  C 

Associating  Zj  with  event  A,  Zj  with  event  C,  i=l,...,k  the  distribution 
theory  proceeds  similarly  to  that  in  Section  2  leading  to  marginal  posterior 
density  estimates  for  p^  and  analogous  to  (4).  Unfortunately  the  restriction 
p.  <  q.  precludes  explicit  calculation  of  complete  conditional  moments  so  that 
there  is  no  analogue  to  (5);  we  use  the  p(r],  s=l,...,v,  to  obtain  posterior 

1 ,5 

*  * 

expectations.  Posterior  density  estimates  for  FA(t  )  and  F^(t  )  along  with 
interpolation  formulas  for  posterior  expectations  can  be  developed  paralleling  the 
discussion  below  (5)  leading  to  (6).  For  the  prior  (11)  distribution  theory  and 

sampling  parallel  that  above  (7).  Posterior  density  estimates  analogous  to  (7) 

*  * 

arise  along  with  posterior  density  estimates  for  FA(t  )  and  for  F^t  ).  The 
use  of  (11)  is  illustrated  with  an  example  in  Section  5 

It  is  worth  noting  that,  regardless  of  prior,  missing  data  can  be  readily 
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handled.  For  instance,  should  X-  be  unavailable  we  need  only  include  it  as  an 
additional  parameter  in  the  model.  Complete  conditional  distributions  for  each 
of  the  other  parameters  when  Xj  is  given  as  well  will  be  exactly  as  before  while 
the  complete  conditional  distribution  of  X^  will  be  Bi(Yj,  Pj/qj). 

5.  A  Second  Numerical  Example 

'Hie  data  in  Table  2  is  taken  from  Maxwell  (1961,  p.  70).  In  his  table 
eacn  of  223  boys  is  classified  into  one  of  five  different  age  groups  and  is  assigned  a 
rating  on  a  four  point  scale  (—1,  0,  1,  2)  for  the  symptom  of  "disturbed  dreams". 
A  rating  of  2  denotes  the  most  severe  suffering  from  bad  dreams  while  a  rating 
of  -1  denotes  no  suffering  at  all.  To  simplify  the  illustration  we  combine  the  two 
intermediate  categories,  ratings  0  and  1,  into  one  to  obtain  Table  2. 

In  the  context  of  Section  4  we  envision  two  underlying  "potency"  curves. 
One  is  associated  with  the  event  A,  "rank  -1  is  assigned",  that  is,  no  suffering 
from  bad  dreams.  The  other  is  associated  with  the  event  C,  "rank  -1,  0  or  1  is 
assigned",  that  is,  no  suffering  or  mild  suffering  from  bad  dreams.  With  these 
definitions,  curves  for  the  incidence  of  A  and  of  C  which  are  monotonically 
increasing  with  age  seem  reasonable.  Obviously  A  c  C  thus  ordering  the  curves. 
Labeling  the  five  age  categories  from  youngest  to  oldest  as  0,  1,  2,  3,  and  4,  for 
category  i,  X^  counts  the  incidence  of  A,  Yj  counts  the  incidence  of  C.  The 
top  part  of  Table  3  converts  Table  2  to  this  notation. 

The  product  Beta  prior,  (11),  is  taken  here.  For  illustration  we  assume 
here  that  all  prior  precisions  M.  are  equal  to  10,  that  FA  is  N(2,4)  and  that 
F  is  N(l,4).  For  each  age  category  prior  incidence  probabilities  associated 

v  ,0 
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with  these  two  curves  appear  in  Table  3  along  with  c^,  and  <5j. 

Table  3  next  lists  the  Bayes  estimates  (posterior  means)  based  upon  r=20 
iterations  and  v=1000  replications.  As  in  Table  1  standard  deviations  and  95% 
equal  tail  interval  estimates  for  each  p^  and  are  also  given.  Lastly  this  table 
gives  the  maximum  likelihood  estimate  for  (p,  q).  The  maximum  likelihood 
estimate  for  p£  badly  violates  the  monotonicity  assumption  but  a  standard  two 
sample  test  reveals  that  it  is  not  significantly  smaller  than  that  for  p^  even  at 
the  0.1  level.  As  analternative  estimate  we  include  the  isotonic  regression  of  the 
maximum  likelihood  estimate  using  the  pooled— adjacent— values— algorithm 
(Robertson,  Wright  and  Dykstra,  1988).  This  estimate  is  still  unsatisfying  since 
strict  monotonicity  might  be  anticipated. 

6.  Summary  and  Discussion 

Nonparametric  bioassay  is  inherently  attractive  in  freeing  us  from 
parametric  specification  of  a  form  for  the  presumed  underlying  potency  curve. 
Similarly,  Bayesian  inference  is  attractive  for  the  bioassay  problem  since  we  often 
have  prior  knowledge  regarding  the  potency  curve,  e.g.,  the  LD50,  the  steepness, 
etc.  While  Bayesian  nonparametric  problems  are  typically  quite  hard,  the  fact 
that  the  bioassay  setting  provides  a  binomial  likelihood  enables  simplification 
particularly  with  regard  to  prior  specification.  Even  so,  previous  attacks  have 
enjoyed  limited  success  in  developing  desired  marginal  posterior  distributions  and 
their  features  due  to  computational  difficulties.  However,  the  Gibbs  sampler 
approach  is  shown  to  be  well-suited  for  the  general  handling  of  such  "likelihood  * 
prior"  forms.  Moreover,  a  potentially  useful  extension  can  be  handled  as  well 
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through  this  sampling-based  approach.  Other  extensions,  for  instance,  to 
multivariate  dichotomous  response  and  to  nominal  polytomous  response  are 
worthy  of  further  investigation. 
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Analysis  of  Data  from  Cox  &  Snell,  1989.  All  n^  =  30. 
Disturbed  Dream  Data  Adapted  from  Maxwell  (1961). 

Analysis  of  Data  From  Maxwell,  1961,  See  Section  5  for  Details. 
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Disturbed  Dream 
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Fig.  1.  Bayes  estimates  of  the  potency  curve. - =  Bayes  estimate  for 

M=l,  •  •  •  =  Bayes  estimate  for  M=100, - =  Prior  curve,  F  , 

while  squares  denote  maximum  likelihood  estimates  at  t=0,l,2,3,4 
Fig.  2.  The  Bayes  estimate  of  F  for  M  =  1  and  its  95%  confidence  band 
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